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Inspired by the work of Kamenev and Kohn, we present a general discussion of the two-terminal 
dc conductance of molecular devices within the framework of Time Dependent Current-Density 
Functional Theory. We derive a formally exact expression for the adiabatic conductance and we 
discuss the dynamical corrections. For junctions made of long molecular chains that can be either 
metallic or insulating, we derive the exact asymptotic behavior of the adiabatic conductance as 
a function of the chain's length. Our results follow from the analytic structure of the bands of 
a periodic molecular chain and a compact expression for the Green's functions. In the case of an 
insulating chain, not only do we obtain the exponentially decaying factors, but also the corresponding 
amplitudes, which depend very sensitively on the electronic properties of the contacts. We illustrate 
the theory by a numerical study of a simple insulating structure connected to two metallic jellium 
leads. 

PACS numbers: 



INTRODUCTION 



The theory of charge transport is often revisited these 
days, half century after the first rigorous quantum formu- 
lation of transport by Kohn and Luttinger,^'^ and after 
the seminal work of Landauer.^ The on going interest is, 
of course, propelled by the recent advances in transport 
measurements at the molecular scale, where all the ele- 
ments that make this problem extremely difficult, from 
the theoretical and experimental point of view, become 
equally important. The challenge comes from the fact 
that we are dealing with out-of-equilibrium, open, inter- 
acting, quantum, many-electron systems. Some recent 
progress in the field includes the derivation of the Lan- 
dauer formulas by Kamenev and Kohn using a linear re- 
sponse (as opposed to a scattering) approach,^ the dis- 
covery of a previously unknown dynamical correction to 
the resistance by Sai and co-workers, the development of 
a master equation approach to transport by Burke, Car 
and Gebauer,^ the elucidation of the role of non local ex- 
change and correlation effects on transport by Burke and 
Koentopp^ the recent discussion of the non-equilibrium 
currents by Doyon and Andrei,^ and the exact solution 
for transport through a quantum dot given by Mehta 
and Andrei.^ All these results remind us that quantum 
charge transport is still a work in progress, even at the 
very fundamental level. 

The recent work of Kamenev and Kohn^ introduced 
an alternative approach to the scattering method for the 
calculation of the two- and four-terminal conductance 
of molecular devices. Although developed within the 
Hartree approximation, the Kamenev-Kohn approach 
can be straightforwardly implemented within the frame- 
work of Time Dependent Current-Density Functional 
Theory.^ We complete this step in the first part of the 
paper, where we show that it is possible to obtain a 
formally exact expression for the adiabatic two-terminal 
dc conductance. In addition, we discuss the dynamical 



contributions to the conductance. When the transverse 
currents are neglected, we recover the dynamical correc- 
tion to the adiabatic conductance discussed by Sai and 
coworkers.— We argue that the transverse currents may 
add further dynamical corrections to the conductance. 

The rest of the paper focuses on the adiabatic two- 
terminal dc conductance. We are interested in the trans- 
port characteristics of long atomic/molecular chains, ei- 
ther metallic or insulating, attached to metallic leads. 
Our goal is to derive the exact asymptotic behavior of the 
conductance with the chain's length. This study comple- 
ments the analysis of Kamenev and Kohn, who focused 
on short junctions. The technique that we use in this 
paper was previously employed in Ref. [ll| to study the 
asymptotic behavior of various perturbations of the elec- 
tron density in metals and insulators, and it is based on 
the analytic structure of the bands and a compact ex- 
pression for the Green's function, as discussed in Ref. [l2|. 
In the present paper, we add a new interesting object, a 
generalized Wronskian with some special properties that 
is very useful in evaluating the conductivity tensor of our 
systems. 

For insulating chains, we obtain the usual exponential 
decay behavior of the conductance with the length of the 
chains. Depending on the structure of the valence and 
conduction bands of the chain, we find that the exponen- 
tial behavior can be associated to more than one relevant 
exponentially decaying term. Carbon nanotubes are typ- 
ical examples of a system in which several exponentially 
decaying terms are important. We find that the asymp- 
totic behavior of the conductance with the chain's length 
is not just exponentially decaying, but may contain mod- 
ulating factors. The exponential decay constants and the 
period of the oscillations can be predicted from complex 
band structure calculations. In addition, we also obtain 
the exact expression for the amplitudes of the exponen- 
tially decaying terms. 

Our work makes a rigorous connection between com- 
plex band structure and tunneling. This connection was 



made in Ref. |lj and, since then, it was further discussed 
in several theoretical and computational studies, some 
of which are mentioned later in our paper. However, 
although the existence of such connection is now quite 
obvious, rigorous and explicit expressions for the ampli- 
tudes of the different exponentially decaying terms were 
missing. Ref. 14 proposes, for example, that the ampli- 
tudes are directly proportional to the number of states 
at the complex k vector. This is only partially true. In 
Ref. [isl, the amplitude is simply set to unity without 
supporting arguments. The explicit expressions for the 
amplitudes, given in the present paper, show that they 
are determined by the overlap integral between the spec- 
tral kernel, the complex k Bloch functions of the chain 
and a suitably defined potential. The amplitudes are de- 
termined by the electronic properties in the immediate 
vicinity of the contacts. The expressions are intuitive 
and simple enough to allow for estimates, without the 
need for costly numerical calculations. 

We also derive explicit expressions for the adi- 
abatic conductance of metallic chains, which show 
the well known oscillatory behavior with the chain's 
length J^iiiii^ii^i^o We discuss the phase and the wave- 
length of the oscillation in the case of monovalent atom 
chains in the light of our analytic results. 

In the last section of the paper, we present a numeri- 
cal implementation of these ideas to a simple insulating 
chain. We illustrate all the elements entering our expres- 
sions for the conductance. For example, we compute and 
display the Riemann surface of the bands. 



II. CONDUCTANCE 

We consider a charge transport experiment involving 
a molecular chain attached to metallic leads, as shown in 
Fig. [TJ Following Kamenev and Kohn,^ we derive a gen- 
eral and formally exact expression for the two-terminal 
dc conductance of such system. 

We start from the Vignale-Kohn linear response 



equation^ 



r'(r,r^c^)Af(r^c^)dr' 



(1) 



which gives the expectation value of the current when the 
system is coupled to a time oscillating electromagnetic 
field. As in Ref. [4, we obtain the dc regime by letting co 
go to zero. 

In Eq. ([1]), x^^{r^r']u;) is the equilibrium Kohn-Sham 
current-current correlation tensor and Af^(r,a;) is an 
effective vector potential. Given the particular gauge 
choice, we can define an effective electric field as Ef^ = 
dtAf^^^ in which case the linear response equation be- 
comes 



j{r,u;)= I a-(r,r';u;)Ef(r',u;)dr'. 



(2) 
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FIG. 1: The geometry of the chain+leads structure. 



A simple and explicit expression for (r', oj) is given in 
Ref. M. 



e« . 



(3) 



where (j)^^"^ is the adiabatic contribution, i.e. the 
linearized static Hartree-exchange-correlation potential, 
and E^^" is the dynamical part of Ef^, given by 



E' 



dyn 



eno 



(4) 



Here, ( is the viscoelastic stress tensor. 

^Ks ^ Kohn-Sham conductivity tensor. 

In the limit ^ 0, the conductivity tensor reduces to 

cT^|(r,r') = ^Tr {j„(r) jf,{r') G^} , (5) 

where = e^? ± Gf^ is the Green's function of the 
equilibrium Kohn-Sham system. 



Gf = (e - ifKs)-\ 



(6) 



and j is the current operator. A convenient expression 
for a is 

<T-|(r,r') = ^Gg(r,r') dl ^ G-(r',r), (7) 

< > > < 

where we used the shorthand da = da— da- An im- 
portant property of the Kohn-Sham conductivity at zero 
frequency is 



^a«<T^|(r,r') = ^Va^(r,r') = 0, 



(8) 



which follows either from the continuity equation applied 
to Eq. ^ or directly from Eq. ((Tj). 

Because of the spatial confinement, the Kohn-Sham 
conductivity tensor goes rapidly to zero as one moves 
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laterally away from the chain- leads structure. We can 
then consider the chain+leads system inside a tube that 
is large enough that the conductivity tensor is practically 
zero at the tube surface and beyond. The net current 
flowing through the molecular chain is given by 



dr' CT'^^(r,r')Ef(r'), 



(9) 



where H is an arbitrary transversal section and the in- 
tegral over is taken only inside the tube. We break 
the current in Eq. (P as / = /"^ + where /"^ is the 
current resulting only from the adiabatic part of the ef- 
fective electric field, Ef = V(0r + ^i''"') = V^f , and 
I"^^"" is the current resulting from the dynamical part of 
Ef. 

To get a clean expression for the conductance, one 
needs to pull out of the integral the physical electric po- 
tential drop between points at z = ±00. This would 
be straightforward if one could make the simplifying as- 
sumption that the effective electric field is uniform in the 
lateral direction. This is, however, a gross approximation 
for the structure in Fig.[Tl The difficulty is not present in 
Ref. 22, which considers a different linear response equa- 
tion, involving the full many-body conductivity tensor 
and the external field. Since one has control on the ex- 
ternal field, it can be considered uniform in the lateral 
direction, greatly simplifying the issue. 

Let us first define a conductance for I"""^ and then com- 
ment on the xc contributions. I"""^ is given by: 

= J dsj dr' (7^'(r,rOV(^f (r ). (10) 

Due to Eq. ([8]), 

a-^(r,rOVVf (rO = Va-^(r, rO(/>i (r^, (11) 

which allows us to transform the volume integral over 
in Eq. (p!Q|) into a surface integral. First, we consider this 
integral over a finite volume, between the surfaces of 
Fig. [H and then take the infinite volume limit by moving 
the surfaces at z = ±oc. An integration by parts in 
Eq. ([9]) gives 

I dsJf - [ I dS'f, <7^^(r,r'X(r'). (12) 

Next we deform the sections I]± into surfaces of constant 
potential. 



(13) 



This is possible because Ti± are arbitrary sections, which 
are used here only to take the infinite volume limit. Then 
it follows that 



P' = r^ f dSa I dS'p cjl}{y,v') 
~4>- I dSa [ dS'p all(v,v'). 



(14) 



After pulling the potential out, the integrals become in- 
dependent of the surfaces, due to Eq. (j8|). Therefore, if 
we deform the sections so that S lies on the xy plane at 
some arbitrary z and Tj± on the xy plane at some ar- 
bitrary z' and take the infinite volume limit, we finally 
obtain 



(15) 



It was argued in Ref. 6 that, within the commonly used 
density functionals, the xc contribution to A(/)^ is iden- 
tically zero. The argument follows from the observation 
that 



Sn 



■ni 



2;=+co 



Sn 



■ni 



(16) 



and the fact that the perturbed density ni is localized 
near the junction. Thus we conclude that is in 

fact the physical electric potential drop A(j)oo and conse- 
quently I"""^ = ^oA(/)cx), where go is a constant depending 
solely on the equilibrium properties of the the system: 



go 



(17) 



This also shows that the conductance derived from the 
adiabatic approximation of the time dependent xc poten- 
tial is exactly the of Eq (pT|) . Thereafter, we call g^ 
the adiabatic conductance. Note that the expression for 
go is formally identical to the one derived by Kamenev 
and Kohn within the Hartree approximation,^ or the 
one derived by Baranger and Stone for non-interacting 
electrons. 

For the net current, we can conclude at this point that 
/ = (^0 + I"^^^ / ^4^00)^4^001 so we can write the total 
conductance g as the formal sum g = go -\- ^^y", where 
gdyn ^ /^y"/A(/)oo. If we neglect the transversal part of 
E^^"" and write it as the gradient of a dynamic potential 
4^1^''^ we can follow the same steps as we did for I"""^ and 
prove that I"^^"" = goA"^"". In other words. 



(18) 



which is precisely one of the key equations (see Eq. 12) 
in Ref. 10. We are then led to conclude that there is an 
additional dynamical correction to the one discussed in 
this reference, correction that comes from the transverse 
part of El^^ 

To evaluate one needs to solve the self-consistent 
equation ([2]) for the current and then compute I"^^"". It 
seems unlikely that one could carry analytic work on g"^^"^ 
beyond that of Ref. 10. We succeed, however, in deriving 
analytic expressions for go in several cases that comple- 
ment the work of Kamenev and Kohni^ 



III. STRICTLY ONE DIMENSIONAL SYSTEMS 

In order to formulate a strategy for calculating the adi- 
abatic conductance, it is useful to start with the simple 
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case of a strictly one dimensional system. In this case, 
the expression for go simplifies to: 

9o = j^G^+{x,x')Kd^'G^-{x,x'). (19) 

The Green's function can always be expanded using 
eigenvectors and eigenvalues. However, the large num- 
ber of terms generated by such expansions are hard to 
control analytically. Our aim is to express go in terms of 
a small number of parameters with straightforward and 
intuitive physical meaning. For this, we use a compact 
expression for the Green's function. 



Ge{x, x') 



(20) 



where x</x> = min / max(x, x') and ^l)f^^{x) are the 
solutions of the Schroedinger equation at energy e, satis- 
fying the boundary condition either to the right or to the 
left. For our infinite system, these boundary conditions 
are simply iljf^^{x) ^ as x ^ Too, respectively. We 
always evaluate the Green's function at an energy e out- 
side the allowed energy spectrum or, at most, take the 
limit of the Green's function for e approaching the energy 
spectrum. It follows from the standard theory of second 
order differential equations in ID that the solutions ifj^ 
are uniquely defined by these boundary conditions. 

One clarification is needed here. When we talk about a 
solution of the Schroedinger equation we do not mean an 
eigenfunction. An eigenfunction is a solution satisfying 
simultaneously the boundary conditions to the left and 
to the right. While the Schroedinger equation has solu- 
tions at any energy, eigenfunctions exist only for certain 
energies which define the energy spectrum. 

W{il)^ (j)) denotes the Wronskian of and ^, i.e. 



(21) 



If il) and (j) are two solutions of the Schroedinger equation 
at the same energy, the right hand side of Eq. ([2T]) does 
not depend on x. Using the Wronskian, we can rewrite 
the adiabatic conductance as: 



^W{i,>i,>_)W{,l,<,i,<) 



(22) 



In the next subsections we evaluate this expression for 
several cases of interest. 



A. Small junction and clean, long metallic leads 

For simplicity we assume that the junction perturbs 
the periodic potential of the clean leads only inside a 
finite interval [0,L]. We denote by h the lattice con- 
stant and by il)k{x) the Bloch functions of the leads. The 
wave number k is complex in general, except when e is 
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FIG. 2: (a and a') The position of in the complex energy 
plane and (b and b') the corresponding kp in the complex k- 
plane for metals (left) and insulators (right). The figure also 
shows the two branch points (red dots), one at /c = and one 
at the zone boundary, and their corresponding energies. 



on the real axis, inside the allowed energy bands. As a 
convention, k will always be chosen with positive imagi- 
nary part {—k will then always have negative imaginary 
part). In ID and only in ID, for each energy e (not nec- 
essarily real), there are two and only two independent 
Bloch solutions, V^/e(x) and i/j-kix). The Bloch solutions 
are normalized according to 



{x)^p-kix)dx = 1, 



and their phase is fixed by requiring 



(23) 



(24) 



One can use the fundamental property of the Bloch so- 
lutions. 



iljk{x^b) 



ikb 



Mx), 



(25) 



to test their asymptotic behavior. One finds that ilJk{x) 
decays to zero when x ^ oo while it diverges when 
X ^ —oo and ip-k{x) decays to zero when x ^ —oo 
while it diverges when x ^ oo. There is an exception 
to this behavior, namely when the energy is in the al- 
lowed energy bands. In this case, the Bloch functions 
behave like waves. Instead of using band indices, it is 
more convenient to think that the wavenumber k lies on 
a Riemann surface 

When the clean leads are perturbed by the junction, 
we have: 



i^fix) - 



and 



ijf{x) 



ipk{x)+7i {k)ip-k{x), X < -L 

T{k)^Pkix), x>0 

T{k)xl;_k{x), x<-L 

tp-k{x)+n+{k)tpk{x), x>0 



(26) 



(27) 



5 



where k is the unique wavenumber with lm(/c)>0 such 
that e/e=e. T{k) and 11^ {k) in Eos. [2611271 are the trans- 
mission and reflection coefficients of the junction. 

To compute the conductance (Eq.[22]) we need the solu- 
tions [26l|27l at Cp. eF is located inside an allowed band of 
the leads and is immediately above/below the allowed 
band, as shown in panel (a) of Fig. [2l The correspond- 
ing wavenumbers k^ are shown in panel (b) of the same 
figure. In the limit J ^ 0, k^ = —kp = /cf, we obtain: 

^±kF{x) ^1Z~{±kF)^^kF{x), x<-L 
T{±kF)^p±kF{x)^ X > 



and 



(28) 



X < -L 



ip^kpix) +Tl+{±kF)tp±kF{x), x>0 



(29) 

By taking x in the right lead, and noticing that 
T{—k) = T{ky, we can easily compute a ffist set of 
Wronskians: 



W{i^>+,^>.) = \T{kF)\^Wo 

^ F ^ F 

W{^<,i>>) = -T{kF)Wo. 



(30) 



By taking x in the left lead, we can calculate the remain- 
ing Wronskians appearing in Eq. ([22|) : 



-\T{kF)\^Wo 
W{4><.,i;>)=T{-kF)Wo, 



^F 

^F ^F 



(31) 



where Wo = ^(V^/cf^ V^-feF)- After straightforward can- 
cellations, Eq. ([22|) leads to the classic result)^ 



5o = -|T(MI^ 

TT 



B. Long chains 



(32) 



Here we consider the opposite case of a clean, long 
molecular or atomic chain attached to arbitrary metallic 
leads. This case is especially important for understand- 
ing the experiments on self-assembled monolayers (SAM) 
grown on a clean substrate. In these experiments the con- 
ductance is probed by an STM tip or a Hg droplet placed 
above the ^AM?^^^^ In this setup, the left and right 
"leads" are very different and the picture in which the 
junction acts as a scattering center for the states of the 
leads is no longer appropriate. 

Deferring discussion of a realistic effective potential to 
the following sections, we assume here a long but finite 
molecular chain attached to metallic leads that is de- 
scribed by an effective potential shown, qualitatively, in 
Fig.[3l The main simplification is that, inside the interval 
[— L/2,L/2], the effective potential is strictly periodic. 



^eff 












X 



-U2 



L/2 



FIG. 3: Schematic of the effective potential of the ID chain. 



By extending this periodic potential to ^oo, one can 
calculate the corresponding energy bands. If the Fermi 
energy of the system, i.e of the infinite leads plus the 
finite molecular chain, falls inside an allowed band we 
call the corresponding molecular chain metallic. We call 
it insulating otherwise. 

Exploiting the fact that Eq. (p!9|) is independent of x 
and x', we conveniently fix x and x' deep inside the chain. 
Then we can use the Bloch solutions of the molecular 
chain in a similar way to what we did in the previous 
subsection with the Bloch solutions of the leads. Nor- 
malization and phase of the Bloch solutions are fixed in 
the same way. 

Inside the [—L/2, L/2] interval, the exact solutions of 
the Schroedinger equation are 



(33) 



where k is again the unique wave number in the upper 
complex semi-plane such that Ck = e. 7^l/r(^) are the 
reflection coefficients of the left/right contacts. The so- 
lutions formally look the same for metallic and insulating 
chains, but the meaning of k^ and the behavior of the 
reflection coefficients at these wavenumbers are different 
in the two cases. 

Metallic chains. In this case ej? is inside an allowed 
band and we can follow the notation in panels (a) and 
(b) of Fig. [21 Taking the limits e ^ and S ^ we 
have: 



(x) = ^l)±kF {x) + 1Z^{±kF)i^^kF (x) 

^F 
^F 

leading to: 

w{i;>,i;>) = [i-\n^{kF)?]Wo 

^F ^F 

w{i;<,i;<) = -[1 - \nukF)\^]Wo 



(34) 



(35) 



and 



w{i^<+,i;>+) = -[i-n,{kF)n^{kF)]Wo 

^F ^F 

w{^<_,^>_) = [i-nu-kF)-R^{-kF)]Wo. 



(36) 
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Since 7^l/r(— ^f) = '^l/r(^f)*, we obtain the fohowing 
simple expression for the adiabatic conductance: 



9q 



1 [l-|7^,(fe^)p][l-|7^^(fe^) 



(37) 



The reflection coefficients contain a phase factor that de- 
pends on the length of the chain and is related to our 
choice for fixing the phases of the Bloch functions. In 
Eq. (|37|) we have assumed the phase of the Bloch func- 
tions to be zero at x = 0. If we move the origin at the 
left/right contacts, the refiection coefficients become in- 
dependent of L. This requires a rescaling:^- 



(38) 



The expression for the adiabatic conductance becomes 



9o{L) = 



1 [i-|7e,(fcf)|^][i-|7e^(fcf)|^] 



(39) 



Eq. ([39|) gives the exact behavior of go as a function of 
L. Since this expression applies to arbitrary periodic po- 
tentials and arbitrary lead potentials, we conclude that 
the behavior shown in Eq. (|39|) is universal for ID metal- 
lic chains. The refiection coefficients are numbers be- 
tween and 1. For poor/good contacts, the magnitude 
of the refiection coefficients is close to 1/0, respectively. 

Resonant transport can also be easily understood from 
Eq. (|39|) . Indeed, ^o(^) is large whenever the denomi- 
nator is small, which may happen at specific values of 
kp' Eq. (|37j) provides not only these values, but also the 
width of the resonance. Indeed, if we look at go as a 
function of kp^ we can see that it has poles whenever kp 
is equal to a solution of the following equation in k: 



e''''^n^{k)n^{k) = 1. 



(40) 



This equation has solutions only in the lower complex 
semi-plane, which can be obtained by analytical contin- 
uation of the refiection coefficients in the lower complex 
semi-plane. The real part of these solutions gives the 
resonant values oi kp and the imaginary part gives the 
width of the resonances. We notice that Eq. (|^Q|) deter- 
mines also the resonances of the Green's function. In 
other words, the positions and the widths of the trans- 
port resonances are exactly the same as those appearing 
in the density of states projected on the chain. To actu- 
ally observe resonant transport, the width of a resonance 
must be smaller than the separation between adjacent 
resonances. 

A consequence of Eq. ([39|) is that go oscillates with 
the chain length. Such behavior was reported previ- 
ously, both experimentally and theoretically^^iiiii^ii^i^S 
According to Eq. ([39|) , the wavelength of the oscillation is 
half the Fermi wavelength of the chain. Chains of mono- 
valent atoms, for example, are half filled and go is pre- 
dicted to be different when the chain contains an odd or 



an even number of atoms: 

i[i-|7e,(MI'][i-I^R(MI'] 



50 (i) = 



TT |l-e-^7^L(A:F)7^H(A:F)|2 



(41) 



This is precisely what has been observed. Interestingly, 
several independent numerical simulations showed that, 
for alkali- metal chains ^li^i^^ g^ is larger for an odd num- 
ber of atoms while for noble-metal chains the opposite 
occurs. ^^^^^ The behavior for noble- metal chains is still 
under debate since existing experiments do not seem to 
confirm the prediction, and a new study finds that al- 
ternative scenarios can happen*^ 

On the basis of expression (|39|) . we can easily under- 
stand these phenomena. In the case of alkali chains, the 
projected density of states on the chain (PDOS) displays 
sharp resonances. The PDOS corresponding to each 
resonance integrates to 2. Thus, for an odd number of 
atoms, the Fermi level sits on top of a resonance while 
for an even number of atoms the Fermi level sits between 
two resonances. According to our discussion of resonant 
transport, when the Fermi level sits on top of the reson- 
cance, the denominator in Eq. (jUJ) is small, leading to a 
large value of go- In fact, if the left and right refiection 
coefficients are the same, go takes the maximum allowed 
value of I/tt (a.u.), no matter how bad the contacts are. 
Now assume that we add one more atom to the chain 
so that we have an even number of atoms. The Fermi 
level moves between two resonances. The refiection coef- 
ficients have slow energy dependence and do not change 
much whereas the phase factor in front of them changes 
sign. The conclusion is that now the denominator of 
Eq. (|^T]) takes a maximum value, leading to a minimum 
value of go- 

Chains of noble- metal atoms are more jellium like*^ 
For good contacts, we expect the refiection coefficients to 
be small and real, in which case Eq. (^1]) predicts that go 
is larger for an even number of atoms. This is the behav- 
ior found in several numerical simulations. ^^^^^ However, 
the phase of the refiection coefficients may be quite sen- 
sitive to the details of the contacts ^^i^^ This can explain 
the different behavior observed in the experiments. 

Insulating chains. In this case, the Fermi energy falls 
within an energy gap. Compared to the metallic case, 
we have the following differences: /c^ ^ /ci? as ^ 0, 
but /cf is now complex as shown in panel (b') of Fig. [2l 
The refiection coefficients have branch cuts along the red 
lines in Fig.[2l Consequently, 7^l/r(^^) are different from 
'^l/r(^f)-^^ Thus, when taking the limit for (5 ^ in 
Eq. ([33|) . we obtain: 



{x) = ^l)kF {x) + n^{kp)^l)-kF (^) 
V^<±(x) = ilj-kF^x) +7^L(A:^)V^/e^(x). 



(42) 



Then: 



w{^>,^>) = [n^{k-)-n^{k^)]Wo 



(43) 
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and 

^ F ^ F 

This leads to 
50 



-[l-7^4fc+)7^H(A:+)]Wo 

-[i-Ti^{kp)n^{k^)]Wo. 

4 Im[7^,(fc+)]Im[7^H(fc+)] 

|l-7^,(fc+)7^a(fc+)p ■ 



(44) 



(45) 



Rescaling the reflection coefficients as in Eq. (|38|) so that 
they become independent of L, we finahy obtain 



4Im[7^L(/c+)]Im[7^K(/c^ 



+ Mg-2/3L 



^ |l-e-2/3i^7^,(4)7^K(4)p 



(46) 



with P=lm{kF)' The above expression shows that the 
behavior of go as a function of L is universal, as the 
above results apply to arbitrary periodic chains and lead 
potentials. 

In the limit of very long chains, the left and right con- 
tacts decouple and the adiabatic conductance becomes 



goiL) = -Im[7^,(4)]Im[7^a(fc^ 



-2/3L 



(47) 



The reflection coefficients are directly proportional to the 
local density of states p^^J^ at the contact edges and at 
the Fermi energy. Indeed, according to Ref. |29|, 



T \nD //+M dep/dp 



peA^)\x^TL/2^ (48) 



which allows us to rewrite the adiabatic conductance as 

2 



5o(i) = - 

TT 



dep/dp 



-2/3L 



(49) 



In the next section, the above expression will be general- 
ized to linear molecular chains in 3D. 



IV. MOLECULAR CHAINS IN 3D 

Real molecular or atomic chains are not strictly one 
dimensional, even though they are electronic systems 
highly confined in two dimensions. Compared to the 
strictly ID systems discussed in the previous section, the 
fundamental difference is that a 3D chain has an infinite 
number of linearly independent Bloch solutions at any 
given energy, rather than just two. In one dimension we 
have scattering processes only between k and —k while in 
3D chains we have scattering processes among an infinite 
set of wavenumbers. 

In the limit of very long chains, like the one schemat- 
ically shown in Fig. [H we can decompose the effective 
potential of the chain+leads into a perfectly periodic 
piece, Vb, extending from — oo to +oo, and a difference 
AV = leff — Vq. The periodic potential Vq is constructed 
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-L/2 



L/2 



FIG. 4: Schematic of the 3D chain+leads. The green shad- 
ing indicates the region that is periodically repeated in the 
construction of the periodic potential Vq. 



by periodically repeating the effective potential between 
—6/2 and 6/2 {b being the chain's lattice constant) in the 
middle of the chain (see Fig.[4|). For very long chains, the 
periodic potential constructed in this way should con- 
verge to the bulk effective potential of the chain rigidly 
shifted by an amount required to align the Fermi lev- 
els of the chain and the leads. The entire procedure is 
illustrated with a concrete example in Fig. [6l 

We may regard the effective Hamiltonian of the chain 
+ leads as a periodic effective Hamiltonian, 

Ho = -^V^ + Voir), Voir + be,) = Voir), (50) 

strongly perturbed by the potential AV. The Kohn- 
Sham Hamiltonian of the entire system is then 



H = Ho^AVLir)^AVRir), 



(51) 



where we divided AV into left and right parts. We are 
not making here the simplification that AVl^r be zero for 
—L/2 < z < L/2, but we assume that AVl^r do decay 
fast to zero as we move away from the contacts. That 
is exactly the case in the numerical example shown in 
Fig. El 



A. Green's function for the unperturbed chain 

The key result that allows us to make analytical 
progress in evaluating the conductivity tensor, is a com- 
pact expression for G^=ie — Ho)~^. This expression is 
a generalization of Eq. ([2Q|) and was recently derived in 
Ref.^. In the same reference, it was shown that the mul- 
titude of Bloch functions, ^n,/c5 and their corresponding 
energies, en,/c, 



[-^V^ + Vb(r)]V^n,/c(r) = en,k^n,kir) 



(52) 



are in fact different branches of global, multi- valued func- 
tions i/jk and e/c, defined on a suitable Riemann surface. 
Later on, we will explicitly calculate this surface for a 
simple system. The band indices do not have any mean- 
ing for complex values of k, since the bands touch and 
hybridize as we move into the complex /c-plane. Thus, 
we abandon the band index and retain only the k index, 
being understood that k lives on a Rimemann surface. 
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The Green's function of the infinite molecular chain is 
given by>i^ 



G°(r,r')=E 



ip-k, (r<)ipki (r>) 
idkCki 



(53) 



where {ki} is the infinite sequence of wavenumbers on the 
Riemann surface such that e^^ = e and r</r> = r/r' if 
z < and r</r> = r' /r otherwise. The Bloch functions 
are normalized so that 



1, 



(54) 



0<z<b 

and their phase is fixed as before. 

B. Green's function for the chain+leads 

We calculate the Green's function for the chain+leads 
from 



(55) 



e — Cr^ -f- Lx^l^Lx^, 

where the Te matrix is given by 

Te = AV + AVGeAV. (56) 
We can naturally decompose the T matrix as 



(57) 



Taking r and deep inside the chain, we obtain 

{T^>fc. (r) Vfc, (r') + Ti^iP.k. ir)i'-k, (r') (58) 
+Tlii>k, (r)^-fc, (r') + T^l^-fe, (r)^^ . (r') } 
where 



(59) 



These coefficients depend, of course, on the energy e. 

At this point we have obtained the exact dependence 
of Ge on the coordinates r and r', which is essential for 
evaluating the conductivity tensor. However, the expres- 
sion of the adiabatic conductance in Eq. (pT|) . with the 
conductivity tensor given in Eq. ([7j) , is far too complex to 
be approached directly. The next subsection is another 
key step in our calculation, which introduces a general- 
ized Wronskian with some remarkable properties that are 
essential to overcoming the difficulty. 



C. Generalized Wronskian 

One can generalize the ID Wronskian in several ways. 
Here we choose the following definition: 



W{^,(fy = j dr^ ^{r^,z)d,^{r^,z). (60) 

This Wronskian appears when calculating the adiabatic 
conductance in Eq. (pT|) . using the conductivity tensor of 
Eq. ([7j) and the above expression for the Green's function 
of the chain+leads. The calculation of the conductance 
further simplifies by virtue of the following remarkable 
property: 



W{ipki,4>-ki) = -2idkekiSki,kj, 



(61) 



where {/c^} is the sequence of wavenumbers corresponding 
to an arbitrary complex energy e. 

To prove Eq. (|6T]) . we consider two arbitrary Bloch 
functions, 



-^V2 + Fo]V^/c = e/cV^/c 



(62) 



From the two Schroedinger equations we obtain the fol- 
lowing straightforward identity 



Zi<Z<Z2 

= 2{ek-ek') J dr^ljkil^k', 

Zi<Z<Z2 



(63) 



where zi and Z2 are arbitrary. The integrand in the left 
hand side of the above equation can be expressed as the 
divergence of a vector field. Furthermore, an integration 
by parts leads to 



r ^ — ^ 



= 2(e/e-e/e/) / dr^ljk^l^k'' 

Z\<Z<Zi 



(64) 



We can immediately see that, if e/e = e^/, then 

w{:i\)k,^k')z=z2 = w{iiJk,Mz^z,^ (65) 

i.e. the Wronskian of the two Bloch functions is indepen- 
dent of z. Now we choose zi = and ^2 = 6 in Eq. (|64|) . 
Using the fundamental property of the Bloch functions, 
we finally obtain 



[e^ik+k> - 1] J d^^li^k'd z^^k']z=o 
= 2(e/e-e/e/) / dr^pk'^k'' 

0<z<b 



(66) 



Thus, if Ck = Ck' but k-\-k' ^ 0, the Wronskian of the two 
Bloch functions is zero. Furthermore, by taking the limit 
k' ^ —k we obtain exactly the second line of Eq. (|6T]) . 
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D. Adiabatic conductance 

Insulating chains. To compute the adiabatic conduc- 
tance, we insert the Green's function (Eq. [58]) into the 
expression for the conductivity tensor (Eq. [7j) and com- 
pute Eq. (pT|) . This process generates Wronskians that 
need to be carefuhy counted. In addition, we need to be 
careful when going from to Cp. 

Let {ki} be the sequence of wavenumbers correspond- 
ing to ej. Since ep fahs in an energy gap, ah wavenum- 
bers are complex. When passing to e^, the wavenumbers 
remain the same, but the T matrix is different. This is 
because we cross the continuum spectrum of the leads. 
When T is evaluated at we call it T. With this nota- 
tion, after counting the Wronskians, the exact expression 
for the adiabatic conductance is 



90 



1 ^ Til 



TT 



1 ^ Ti^n 



rfiij rpij 



(67) 



rpt'J rplj _ rplj rpl^ 

J- lrJ- rl J- lrJ- r 



TT 



No approximation has been made so far. 

We now show that Eq. (|67|) greatly simplifies in the 
limit L ^ oo. We introduce the following notation: 



(3i lm{ki), (3min min{A}- 
First we demonstrate that: 

Tlk = [1 + o(e-2^-^"^)]TLG^TK. 



(68) 



(69) 



This follows from the following considerations. First, we 
notice that every time that is sandwiched between a 
AVl and a AV^, it becomes of order o(e~^"^^^^). This fol- 
lows from the exact expression of given above. Next, 
we iterate Eqs. (|55|) and ([56|) to generate a formal series 
for Tlr: 



T,^ = AKG^Ay^ + AV,G^,AVG^,AV^ 



(70) 



+AKG'^AVG^AFG^AFk + . . . 
Let us consider, for instance, the third term in Eq. ([7Q|) : 



+AV^G°AV^G°AV^G°AV^ 
+AV^G°AV^G°AV^G°AV^ 
+AV^G°AV^G°AV^G°AV^. 



(71) 



In the first three terms, G^ is sandwiched between a Al/^ 
and a AFr only once, whereas all the three G^ that ap- 
pear in the last term are sandwiched between a AVi^ and 
AV^. As a consequence, the ratio of the fourth term and 
anyone of the first three terms is of order o(e~^^"^^"^). 



By applying this argument to all the terms in Eq. ([7Q|) 
we obtain: 

n 

^Zk=E AKG^..AK AF,...G^AF, 

k=l V ' ^ V ' W^j 



k 



n -\- 1 — k 



plus terms o(e~^^"^^"^) times smaller. 

Next, we consider the expansion of Ti^G^T^^ in powers 
of G^. By applying the same arguments that led us to 
Eq. ([72|) . we find that the n-th term in the exapnsion 
is equal to the right hand side of Eq. ([72|) plus terms 
Q(^^-'2f3minL^ times smaller. This proves Eq. (|69|) . Taking 
the matrix elements of this equation and using the rep- 
resentation of the Green's function G^ given in Eq. ([53|) 
gives: 



Tli = [l^o{e- 



L 



(73) 



Similar conclusions holds for T^l and for the tilde coun- 
terparts. 

As shown below in Eqs. ([77|) and ([79|) . the matrix ele- 
ments Tl^ and T^^ are exponentially small for long chains. 
Thus, if we retain only the leading terms, Eq. (|67|) be- 
comes: 



30 = -Y] 



dk^kidk^kj 



(74) 



This is the exact asymptotic form of qq. The terms that 
we have neglected are o(e~^^"^^"^) times smaller. 

Finally, we show that the matrix elements of T have 
simple and intuitive expressions. For instance, the first 
factor in the numerator on the right hand side of Eq. ([74|) 
is: 

Tli - fli = {^_fcJAK(G4 - G,-)AK|V-fe,). (75) 
This can be expressed in terms of the spectral operator 



Pe 



(76) 



The diagonal elements, pep{x^x)^ of the spectral opera- 
tor give the local density of states. By writing V^/c(r) = 
Uk(v)e^^^ ^ with Uk periodic, Uk(v + he^) = Uk{r), we ob- 
tain 



with 



rpij _ rj^ij _ ^±(ki-\-kj)LQij 



u_k. (r) AK {r)pe, (r, AK {r')u-k. (r^, 



(77) 



(78) 



where r and are measured from the left end of the 
chain, z = —L/2. Similarly 



jiij _ jiij _ ^±(^ki-\-kj)LQij 



(79) 
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with 



(80) 



where r and are measured from the right end of the 
chain, z = L/2. In the hmit L ^ oo, the coeffi- 
cients become independent of L. By inserting Eqs. ([77|) 
and ([79|) into Eq. ([74|) . we obtain the fohowing simple 
asymptotic form for go: 



7..1 * -J 



(81) 



We notice that the integrands in Eqs. ([78|) and ([8Q|) con- 
tain the fohowing terms: 



g-i(/c,^+fe,^')AK(r)AK(rO 
gi(fc,^+fc,^')A14(r)A14(rO, 



(82) 



which are highly localized near the left and the right con- 
tacts, respectively. This shows that the conductance of 
the molecular device is only determined by the properties 
of the chain and of the contacts. 

Strictly speaking, the asymptotic form of ^o(^) is de- 
termined by the wavenumber k such that Im(/c) = /3min- 
However, for complex molecular chains such as carbon 
nanotubes,^^ there may be many wavenumbers with sim- 
ilar imaginary parts, especially when the valence and the 
conduction bands are highly degenerate. 

Metallic chains. As in the insulating case, let us con- 
sider the infinite sequence of wavenumbers {ki} for which 
^ki = ^F- Since the Fermi energy falls within allowed 
bands, some of the ki are real. We can restrict ourselves 
to the latter, because complex wavenumbers lead to ex- 
ponentially small contributions to go in the limit L ^ oo. 

In order to compute the conductivity tensor in Eq. ([7j) 
we need both G + and G - . When going from 



jp to e^. 



ki becomes — fc. Thus: 



T^j ^ f^j = (r^r, dkCk, ^ -dkCk,. (83) 

Following the same path as in the insulating case, we 
obtain: 



90 



1 



2ImT.*: 



(84) 



TT 



dk^kidk^kj 



At zero temperature, the left and right contacts do not 
decouple and all the terms in Eq. ([84j) should be retained, 
even in the limit of very long chains. However, at finite 
temerature, the Fermi energy acquires a small imaginary 
part and the same decoupling mechanism of the insulat- 
ing chains holds. In this regime, Eq. ([84|) simplifies to: 



90 



N 

TT 



y 



(85) 



(a) 
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FIG. 5: Upper panel: The self consistent electron density. 
Middle panel: The self consistent potential. Lower panel: 
The local density of states, integrated over the xy cooridnats. 
The black line shows the position of the Fermi level. The 
numbers represent atomic units. 



where N is the number of real k points such that Ck = ep- 
If the potential Vo(v) is separable in longitudinal and 
transverse coordinates, Vb(r) = Vo{x^y) + Vq (^), N is 
simply the number of transverse modes. ^ In the general 
case, N is the number of eigenstates of the T matrix at 
e = ep' In our approach, N can be extracted from the 
band structure of the chain. For instance, with reference 
to the band structure in Fig. [HI if ej? falls within the 
lowest energy band, N = 1. li ep falls within the third 
band (in order of increasing energy) and crosses the band 
at more than two points, N = 2. 

The matrix elements of T in Eq. ([84|) exhibit oscilla- 
tory behavior with L. Thus the same conclusion that 
we reached for strictly ID case remains valid, namely, go 
should exhibit oscillatory behavior with L. Again, the 
phase of the oscillations will be affected by the resonant 
behavior of the T matrix. 



V. NUMERICAL EXAMPLE 

In this section we provide an example that illustrates 
how the above ideas can be implemented numerically. 
We consider an insulating chain made of quantum dots 
(qdots). These are jellium qdots modeled by a spherically 
symmetric potential well: 



-vo 



1 + exp[(r - ro)/0' 



^0 



O.la.u. 



(86) 



The smoothness of the potential edge is controlled by ^, 
which we take here to be ^ = 1 a.u.. We also added a 
neutralizing background of positive charge, with charge 
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FIG. 7: Left panel: The lowest bands of the periodic potenial. 
Right panel: The Riemann surface corresponding to these 
bands. 



FIG. 6: Upper panel: The effective Kohn-Sham potential of 
the chain + leads and the region that is repeated periodi- 
cally. Middle panel: The periodic potential obtained by re- 
peating the central piece of the upper panel. Lower panel: 
The difference between the effective and periodic potential. 
The numbers represent atomic units. 



density: 



n+(r) = 



no 



1 + exp[(r - ro)/0' 



(87) 



where no has Ts = 2.5 a.u. and tq = 3.15 a.u.. These val- 
ues were chosen so that n+ integrates to 2e. Thus, the 
qdot can accommodate 2 electrons and has a complete 
shell. A linear periodic array of such qdots will be insu- 
lating if the distance b between adjacent qdots is large 
enough. This condition is met by our choice b = 12 a.u.. 

We have considered a chain of seven qdots attached to 
cylindrical jellium leads at both ends. The setup is shown 
in Fig. [SJa). Periodic boundary conditions were adopted 
along z (with a cell parameter of 124 a.u.), while free 
boundary condition were used in the lateral directions. 

In Figs.[5][b) and (c) we display the self-consistent den- 
sity and potential Kff(i*), respectively. Fig. [5](d) shows 
the local density of states integrated over the xy coordi- 
nates. Here one can see the bands of the chain and the 
Fermi level, which is located between the first two bands. 

Fig. [6] illustrates the procedure that we use to construct 
the periodic potential Vb(r). This is obtained by repeat- 
ing infinite times the part of Kff(r) comprised between 
—b/2 and b/2. The resulting periodic potential is shown 
in Fig.[6](b). Fig.[6](c) shows the difference AV = Ven—Vo. 
The potential is cylindrically symmetric, so the Lz quan- 
tum number is conserved. Since the lowest bands corre- 
spond to the Lz = sector, the entire discussion will be 
restricted to this sector. 

The next step is the complex band structure calcula- 
tion and the construction of the Riemann surface for the 
periodic potential Vo{r). The theory of the complex band 
structure for a general periodic potential was discussed in 
Ref. tl2i. Based on this theory we compute the Riemann 



surface of the bands for Vb(r). In Fig. [71 we show the 
four lowest bands and the corresponding section of the 
Riemann surface, which consists of four unit disks. The 
Riemann surface is represented using the natural vari- 
able A = e*^^. The full dots indicate branch points. The 
dashed lines connecting the branch points to the center 
of the disks represent branch cuts, which always occur in 
pairs. In Fig. [Tithe pairs are located on adjacent disks. 
The disks become connected via these branch cuts. 

In order to construct the Riemann surface, we diago- 
nalize the k • p Hamiltonian, 



i(-iS, + fc)2+KeO (88) 



defined in the interval [—6/2,6/2] with periodic bound- 
ary conditions. The k • p Hamiltonian was implemented 
on a real space grid, using a 5 point finite difference rep- 
resentation for the second derivatives. The grid spacing 
was 0.5 a.u.. 

Each panel of Fig. [51 shows the eigenvalues en,/c, n = 
1, . . . , 4 of the k • p Hamiltonian as functions of /c, when 
the imaginary part of k is fixed and the real part is varied 
from —ir/b to ir/b. When k moves on this line in the 
complex /c-plane, the variable A = e*^^ moves on a circle; 
the larger lm{k) the smaller the radius of this circle. We 
report in Fig. [51 a set of panels going from (a) to (j). Each 
panel corresponds to a different value of Im(/c), increasing 
from lm{k) = a.u. (panel (a)) to lm{k) = 0.12 a.u. 
(panel (j)). Each panel shows two diagrams: the one on 
the left represents Re(en,/c) vs Re(/c), and the one on the 
right represents Im(en,fe) vs Re(/c). 

In order to understand the analytic structure of the 
bands, for a given value of k (or A), one may draw a verti- 
cal line like the one in the right panel of Fig. [71 This verti- 
cal line intersects the Riemann surface infinite times. The 
infinite sequence of complex energies en,k^ = '^^ • - • ^oo, 
of the k • p Hamiltonian are in fact the values of the same 
function, e/e, evaluated at these intersection points. If we 
move A on a circle of radius e~^^^^^\ these intersection 
points generate certain contours on the Riemann surface. 
When e/e is evaluated for all A located on these contours. 
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one obtains certain paths in the complex-energy plane. 
The complex bands of Fig. [8] are in fact examples of such 
paths. 

If a contour remains on the same Riemann sheet as 
A completes the circle, e.g. the n-th sheet, then en,k 
moves continuously on a closed path in the complex en- 
ergy plane. In other words, both Re(en,/c) and Im(en,/c) 
return to their initial value upon completing the circle. 
If this does not happen, i.e. the complex en,k does not 
return to its initial value, the contour has intersected a 
branch cut. When A hits a branch point connecting the 
sheets n and m, the complex bands generated by en,k and 
em,k intersect. 

The Riemann surface of Fig. [71 was computed as fol- 
lows. For small Im(A:), all Cn^k return to their initial value 
when Re(/c) hits the right edge of the Brillouin zone. This 
is what we see in the first two panels of Fig. [HI As we 
increase Im(/c), we hit the first branch point (labeled 1 
in Fig. [7j). This branch point connects the third and the 
fourth sheet, thus the real and the imaginary parts of cs^k 
and /c are equal at this particular value of k. This is 
exactly what we see in panel (c). By further increasing 
Im(/c), we hit the second set of branch points, indicated 
by 2 and 2' in Fig. [7l These branch points connect the 
second and the third sheet, thus the real and imaginary 
parts of 62, /c and 63, ^ are equal at these k points. In panel 
(d) we have crossed these branch points. Indeed, by ex- 
amining the real and imaginary parts of 62,^ and 63, /c, we 
see that they no longer remain on the same sheet. For ex- 
ample, 62, /c moves from the second to the third sheet and 
then to the fourth sheet. Next we hit the third branch 
point, and we see this happening in panels (d) and (e). 
At last, we hit the fourth branch point, which connects 
the first and second sheet. We see this happening in pan- 
els (i) and (j). By further increasing lm{k) we do not find 
any additional branch points on the first three Riemann 
sheets. On the fourth sheet, there are additional branch 
points connecting this sheet with upper sheets. 

Finally, we compute the coefficients. In the case 
that we are considering here, the Fermi level falls be- 
tween the first and the second band. The complex band 
calculations provide the corresponding sequence of {ki} 
wavenumbers. We need to consider only the k wavenum- 
ber located on the branch cut connecting the first and 
the second Riemann sheets. Consequently, we can drop 
the ij indices in Eqs. ([TS]) and (|8Ql) . The complex band 
calculations also provide the left/right eigenvectors U-k 
and Uk that enters the definition of the O coefficients. 
The coefficients are then calculated as follows. Rather 
than computing the spectral operator, we obtain directly 
the functions G^i AVl/rV^^/c, by solving for ^ in 

(e±-if)* = AK/a^/'Tfc- (89) 

Eq. (|89l) was solved on a real space grid with the same 
grid spacing as in the case of k-p equations. The metallic 
leads were increased until the supercell reached 250 a.u. 
in length. With this supercell, we could use an imaginary 
part of S = 0.005 Ha for e^. After we calculate ^, we 




FIG. 8: Complex band calculations. Each panel shows the 
real (left) and imaginary (right) parts of en,fe, ri = 1,...,4, 
for different values of Imk. In the panels from (a) to (f) 
Im(A:) was uniformly increased from to 0.06 (in a.u.). In 
the panels from (g) to (j), Im(A;) = 0.1, 0.11, 0.115 and 0.12 
(a.u.), respectively. 



take its scalar product with 27riA]/L/R'0^/c to finally get 
6l/r from Eqs. ([78l) and (|80l). 

To get more insight into the transport properties of 
the chain, we have rigidly moved the Fermi level by a 
biased potential V and calculated ^0 at 6i? + e<I>. To first 
approximation, this will give the non-linear differential 
conductance of the chain. We have repeated the above 
steps for chains containing 2, 3,. . . , 7 qdots and the re- 
sults are shown in Fig. [9l 

VI. DISCUSSION 

The self-consistent potential is not the same for the dif- 
ferent structures considered in Fig. [9l Consequently, the 
periodic potential used in the calculation of go is different 
from case to case. However, due to nearsightedness,^^ in 
the middle of the chain we expect these self-consistent 
potentials to converge to the bulk value, as the chain's 
length gets larger and larger. In our calculations, we 
have seen a substantial change in the periodic potential 
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FIG. 9: The differential conductance for chain+leads struc- 
tures containing 2 (black), 3 (blue), 4 (green), 5 (red), 6 (cyan) 
and 7 (indigo) qdots, as function of the applied bias poten- 
tial. The thick vertical dashed line represents the converged 
position of the Fermi level. The thin dotted lines indicate the 
edges of the valence and conduction bands for the structure 
containing 7 qdots. 



for the systems containing 2, 3 and 4 qdots. However, 
the band energies of the periodic potentials for the struc- 
tures containing 5, 6 and 7 qdots are almost the same, in- 
dicating that the periodic potential is almost converged. 
The edges of the valence (e^) and conduction (cc) bands 
for the periodic potential corresponding to the structure 
with 7 qdots are indicated by thin, vertical dotted lines 
in Fig. [9l The Fermi level, marked by a thick vertical 
dashed line in Fig. [9l was practically the same for all 
systems. 

Each circle in Fig. [9] represents the output of a single 
numerical run. One can see that go varies smoothly with 
the bias potential. This indicates that the super-cell used 
in the calculation is large enough, so that the density of 
states appears smooth when computed with our damping 
parameter 6 = 0.005 Ha. The points in Fig. [9] can be 
connected by a line having the characteristic shape found 
in Ref. [H. There is one important difference, however: 
our line shape is not parabolic as in Ref. 33, but has an 
asymmetry since the conductance minimum occurs closer 
to the valence band. 

go is primarily determined by the variation with the 
bias potential of Im(/c), dk^k and ipk- For very long 
chains, the conductance minimum occurs at the maxi- 
mum value of lm{k). In our example, the conductance 
minimum is achieved when ep -\- is equal to the en- 
ergy of the branch point, i.e. when k becomes the branch 
point itself. This is not a universal feature, however, be- 
cause the energies of the branch points are complex, in 
general, and cannot be equal to cf + e^- For example, 
this is the case with alkyl chains. The asymmetry of the 



curves in Fig. [9] is in agreement with the theoretical re- 
sults of Ref. 29 about the position of the branch point, 
which is reviewed below. 

When the gap is small relative to the width, of the 
valence band, one can use, quite accurately, the effective 
mass or quadratic approximation for the band energies, 
in which case finding the complex k vectors for all en- 
ergies in the gap is trivial. In particular, the energy of 
the branch point is found to be in the middle of the gap. 
The fact that one can use the effective mass approxima- 
tion to calculate the complex k wavenumbers, which give 
the exponential decay with the chain's length of go, im- 
plies that one can also use a simple square barrier model 
to calculate go. In other words, if one plugs the trans- 
mission coefficient for tunneling of a particle of mass m* 
(the effective mass of the electrons in the chain) through 
a square barrier of height ec — ei? — e<l> in the Landauer 
formula, one will find the correct exponential decay be- 
havior of go with the chain's length. 

In the narrow band limit, which is our case, the branch 
point strongly shifts towards the valence band and its 
imaginary part is found to be:^^ 



bw 



(90) 



This shift leads to the asymmetry that can be seen in 
Fig. [9l For example, for the chain+leads structure con- 
taining 7 qdots, Eq. ([9Q|) gives Im(/c)=0.14 a.u., as op- 
posed to a value of about 3 a.u. that is provided by the 
effective mass approximation. The exact value of lm{k) 
is 0.115 a.u.. This is a concrete example of the break- 
down of the effective mass approximation for chains with 
narrow bands. Because of that, the picture of electrons 
tunneling through a square barrier is misleading in this 
regime. 

Note that dk^k is zero at the band edges. This means 
that our asymptotic expression for go diverges as we move 
towards the gap edges. This can be corrected by includ- 
ing multiple-reflection terms in Tl^. Also note that dk^k 
is infinite at the branch point. Apparently, this leads 
to a singularity. However, the Bloch functions, t/jk and 
tp-k, also diverge at the branch point. The divergences 
cancel exactly and go becomes finite and smooth at the 
branch point. 

Besides the exponential decay, the conductance should 
also display oscillatory behavior with the length of the 
chain. This oscillatory behavior comes from the real 
parts of ki in Eq. ([8T]) . In our example, the oscillation is 
absent because k is at the zone boundary. The oscillation 
would be present if the Fermi level were located between 
the second and the third band of Fig. [71 

For all six chain+leads structures, we found that the 
Fermi level practically coincides with the energy of the 
branch point. This was quite surprising, because the 
branch point in our calculations is not in the middle of 
the gap as in Refs. 34-35. This sends us back to Tersoff's 
criterion, which says that the Fermi level is pinned at the 
branch point. The arguments leading to this conclusion 
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are valid for cases when the metal- induced gap states ex- 
tend far into the chain and are based on the assumption 
of local (as opposed to global) charge neutrality. In ID, 
one can use an observation by Rehr and Kohn, which 
says that the gap states below/above the branch point 
originate from the valence/conduction band.^^ Then, in 
order to maintain the local charge neutrality somewhere 
deep in the chain, we must occupy only the gap states 
from the valence band. This pins the Fermi level at the 
branch point. A clear division between gap states is al- 
ways possible in ID because the energy of the branch 
point is real. This is also true in our example but the 
energy of the branch points connecting the valence and 
conduction bands can be complex. Thus, no such clear 
division of the gap states is possible. Nevertheless, Ter- 
soff 's criterion seems to be confirmed by many numerical 
applications, including ours. It will be interesting to see 
if Rehr and Kohn's observation can be generalized from 
ID to real molecular chains and if a more general expla- 
nation of Tersoff 's criterion can be found. 



VII. CONCLUSIONS 

In conclusion, we derived a formally exact expression 
for the two-terminal dc conductance within Time Depen- 
dent Current-Density Functional Theory. This general 
result may provide a useful starting point to go beyond 
the adiabatic approximation in transport. 

Using a compact expression for the Green's function, 
previously reported in Ref. [l2|, and a generalized Wron- 
skian, reported here for the first time, we derived explicit 
analytic expressions for the adiabatic part of the conduc- 
tance. For insulating chains, we found that the adiabatic 
conductance is proportional to the overlap between the 
spectral kernel pe^(r, r^), the complex k Bloch functions 
of the chain and the potential perturbation at the con- 
tacts AVl/r- We also found that the conductance has 
an exponential decay behavior with the chain's length. 



This behavior may be modulated by an oscillatory fac- 
tor. Both the exponential decay constant and the period 
of the oscillations can be predicted from complex band 
calculations. For metallic chains, we rediscovered the os- 
cillatory behavior of the conductance with the chain's 
length. In the light of the new analytic expressions, we 
discussed the phase differences in these oscillations that 
were observed in alkali and noble metal chains. 

The results of this paper may open the possibility for 
quantitative theoretical predictions on tunneling trans- 
port through extremely long molecular chains that can- 
not be treated by current ab-initio approaches. The re- 
sults can be extended to finite temperatures and thus 
may provide a quantitative understanding of the tunnel- 
ing transport experiments involving alkyl chains grown 
on a silicon surface, in both thermionic and tunneling 
regimes. The analysis can be also generalized to the spin 
dependent case, useful to understanding the tunneling 
magneto-resistance. We think that even in this case it 
may be possible to obtain quantitative predictions sim- 
ilar to those of Ref. i38|, without going through costly 
ab-initio calculations. 
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